{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Asymptotic Bode diagrams\n",
    "\n",
    "One of the big advantages of Bode diagrams is that they are very easy to sketch out by hand (or, equivalently, to visualise mentally)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = numpy.logspace(-2, 2, 1000)\n",
    "s = 1j*omega"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Systems with real poles\n",
    "\n",
    "Let's study the bode diagrams of systems of the form\n",
    "\n",
    "$$\\frac{K}{(\\tau s + 1)^n}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def annotated_bode(ax_gain, ax_phase, G, K, tau, order):\n",
    "    high_freq_asymptote = K/(tau*omega)**order\n",
    "    \n",
    "    # Gain part\n",
    "    ax_gain.loglog(omega, numpy.abs(G))\n",
    "    ax_gain.axhline(K, color='grey')  # Rule 1\n",
    "    ax_gain.loglog(omega, high_freq_asymptote, color='grey') # Rule 2\n",
    "    ax_gain.axvline(1/tau, color='grey')  # Rule 2\n",
    "    ax_gain.set_ylim([1e-2, 1e+1])\n",
    "    ax_gain.set_ylabel('|G|')\n",
    "\n",
    "    # Phase part\n",
    "    ax_phase.axhline(0, color='grey')  # Rule 3\n",
    "    ax_phase.semilogx(omega, numpy.unwrap(numpy.angle(G)))  \n",
    "    ax_phase.axhline(-numpy.pi/2*order, color='grey')  # Rule 4\n",
    "    ax_phase.axvline(1/tau, color='grey')  # Rule 5\n",
    "    ax_phase.set_ylim([-3*numpy.pi/2, 2*numpy.pi/2])\n",
    "    ax_phase.set_ylabel(r'$\\angle G$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotresponse(order=1, tau=1, K=1):\n",
    "    G = K/(tau*s + 1)**order\n",
    "    \n",
    "    fig, [ax_gain, ax_phase] = plt.subplots(2, 1)\n",
    "    annotated_bode(ax_gain, ax_phase, G, K, tau, order)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0a312237c6d1476d9824affe955b23cc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='order', max=3, min=-2), FloatSlider(value=1.0, descripti…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(plotresponse, order=(-2, 3), tau=(0.1, 10), K=(-1., 2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that we can construct a reasonable approximation by knowing a couple of things\n",
    "\n",
    "1. The gain ($K$) of the system defines the low frequency asymptote of the gain graph\n",
    "2. The high frequency asymptote of the gain is $\\frac{K}{(\\omega\\tau)^n}$. Effectively, on a loglog scale, this means we have -n/decade slope above frequencies of around $1/\\tau$\n",
    "3. The low frequency phase asymptote is 0\n",
    "4. The high frequency phase asymptote is $-n\\pi/2$\n",
    "5. The phase curve has an inflection at $1/\\tau$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Systems with complex poles\n",
    "Systems with complex poles show uniique frequency response behaviour. We will focus on the second order system shown below:\n",
    "\n",
    "$$ G = \\frac{K}{\\tau^2 s^2  + 2\\tau\\zeta s + 1}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotresponse(K=1, tau=1, zeta=1):\n",
    "    plt.figure(figsize=(15, 5))\n",
    "    order = 2\n",
    "    G = K/(tau**2*s**2 + 2*tau*zeta*s + 1)\n",
    "    \n",
    "    ax_gain = plt.subplot2grid((2, 2), (0, 0))\n",
    "    ax_phase = plt.subplot2grid((2, 2), (1, 0))\n",
    "    ax_complex = plt.subplot2grid((2, 2), (0, 1), rowspan=2)\n",
    "\n",
    "    annotated_bode(ax_gain, ax_phase, G, K, tau, order)\n",
    "    \n",
    "    # poles\n",
    "    poles = numpy.roots([tau**2, 2*tau*zeta, 1])\n",
    "    ax_complex.scatter(poles.real, poles.imag)\n",
    "    ax_complex.axhline(0)\n",
    "    ax_complex.axvline(0)\n",
    "    ax_complex.axis([-2, 2, -2, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0d52c9d411ea4ec99ff02befcc62ee71",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='K', max=2.0, min=0.1), FloatSlider(value=1.0, descri…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plotresponse(K=1, tau=1, zeta=1)>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interact(plotresponse, K=(0.1, 2), tau=(0.1, 2), zeta=(0., 1.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the rules from before still hold, except that we start seeing the so-called \"harmonic nose\" emerge when $\\zeta<\\sqrt{2}/2\\approx{0.7}$. The maximum of the nose occurs at the resonant frequency of\n",
    "\n",
    "$$\\omega_r = \\frac{\\sqrt{1 - 2\\zeta^2}}{\\tau}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dead time\n",
    "\n",
    "The effect of dead time is to increase the phase lag indefinitely as a function of frequency. Delay has no effect on the gain of a system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = 1\n",
    "G = numpy.exp(-D*s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x117a944e0>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHPZJREFUeJzt3XmYVHed7/H3t6p636G7aegFusPSAkIIHSBK9mjImIgacwV14jhRHCfqTLxz75jrvc7i9bkzk5mrN65DTB6XOMYYjYk6ionRxGgIAYIJSyANBGjWZm26obeq3/2jqknTNIHq6lOnqs/n9Tz19Dm/s3392ekPv7OVOecQEZHgCvldgIiI+EtBICIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAJOQSAiEnARvwu4GJWVlW7KlCl+lyEiklXWrVt32DlXdaH1siIIpkyZwtq1a/0uQ0Qkq5jZrotZT6eGREQCTkEgIhJwCgIRkYBTEIiIBJxvQWBmS8xsq5m1mtln/KpDRCTofAkCMwsDXwVuAmYCy81sph+1iIgEnV+3jy4AWp1zOwDM7CFgKbB5NA/S0x/l6a3tJI4R/5lYlpjFDGyg9ewfw24zsK4NWff1be2cZWf2M6jdhmw43Dbn1Dpk3TdaNlz70P2RqDVkRmjgZ2jQtJ09bSEID0wn2sOJ9W1wUSKSVfwKglpgz6D5NmDh4BXMbAWwAqChoWFEB+ns7mfFd9eNsERJVsggHDLMhgmT0LnBEjIIh42cUIhI2AiHQuSEjUjIiITj0+FQiJyQEQnH2yIhIzKwXjg+PXj9nHCIvEjikxMmPydEXiScaBs0nxMiP/FzYFleJEQopECT4PErCIb7r82dNePcSmAlQEtLixtm/QsqLcjhZ59cPGifAwdyZ+bdmWXurCLcmSMOt+75txnY95AfZx3bnbPs7P3gzq7x7O1f3+aNaj1fXcPux0E05oglpmPOEXWOmIsvj8Xi0zHnEp/E9FntJOYHbecc0RiJ/cbbowPTsfh0NOboi8bojzr6Y47+WHy6Lxqjpy9GXyxKf2J5XyxGNObOLO9PbHumLRYb1BcjkxcJUZwXoTAvTFFuJDEdoTgvTOHAfG6YorwIRblhivNzKCvIobwwh/KCHMoK4/N5kXBqhYikkV9B0AbUD5qvA/aN9kFywiFm15aN9m4lg/VHY/QmQqS7P0pPX4ye/hg9/VG6++I/B9q6+6Jnlg3Mn+6N0tXbT1dPlK6efrp6+zlxuo/9x08n5uPt/bE3TpyCnDDlha+HRFlBDuOK8qgqSXyK4z+rS/KoLM6jIFfBIf7xKwheAKaZWSOwF1gGvN+nWmQMiYRDRMIhCnO9O4Zzjp7+GKd6o3R2x4Pi+Olejp/q4/jpPk6c6o23nZnvY0d7F+t2HeNIV++wo5aSvMiZkKitKKCuvIDaigJqywupqyhgYnm+RhniGV+CwDnXb2afAFYBYeAB59wmP2oRSZaZkZ8TJj8nzLii5BKnPxrjaFcvh0720N7ZQ/vJQZ/OHg6e6Oa57Uc40NF9TmBUl+RRP66QpsoimqqKaaoq4pKqIhrGFZEb0SNBMnLmUj2pmgYtLS1OL52TIOmLxjhwops9x06x99hp9h4/zd5jp9l19BQ72rs43NlzZt1wyKivKGBqdTEzJ5bypomlzJxUSn1FoS5+B5yZrXPOtVxovax4+6hI0OSEQ9SPK6R+XOGwyzu646ebdrR3xn8e7mTbwU6eeuUQA5cvivMiNNeUMHNSKfMayrmsoYKGcYW61VfOoRGByBjS3Rdl64GTbNnfweb9HfGf+zro6o0CUFmcy7yGCi5rqKBlSgVz68p1WmkM04hAJIDyc8LMrS9nbn35mbZozLHt4EnW7z7Gul3HeHH3cZ7YfBCAwtwwCxvH8daplSyeVsmMCSUaMQSQRgQiAXSks4cXXjvGH7Yf5tnWw+xo7wLiI4arp1fztpkTuGp6JYW5+rdiNrvYEYGCQETYd/w0v2+Nh8Jvt7Zz4nQfeZEQV06r5O0za3j7rAmUe3lPrnhCQSAiI9IXjfHCzqP8avNBnth8kL3HT5MTNq6ZUc2759VyXXM1+Tl6piEbKAhEJGXOOTbt6+CxDXt5bMM+Dp3soSQvwjvmTOT9CxuYU1d+4Z2IbxQEIjKqojHHc9uP8OiLe/nFxv2c6o0yp66MDy6azC1zJuk1GRlIQSAinuno7uPR9Xt5cPUuXj3USWl+hOULG7hjcSPVJfl+lycJCgIR8Zxzjud3HuW7z+3iFxv3EwmHeO/8Oj52VROTxxf5XV7g6TkCEfGcmbGoaTyLmsbz2uEu/v2ZHTyyto2H1uzmnXMncdfbpisQsoBGBCIyqg51dHP/szv59nOv0R91LF/QwCevn6pTRj7QqSER8dWhjm7ufepVHlqzh5xwiI9e2cjHr5mqi8ppdLFBoJeMiIgnqkvz+d/vejNPfvpqbpg5gXufauWG//s0v9x4gGz4B2iQKAhExFNTKov48vJ5PPyxKyjJj/AXD67j9gfWsOtIl9+lSYKCQETSYkHjOH72ycX83S0z2bD7OEu+9Du+9fudxC7wtZ/iPQWBiKRNJBziw29t5IlPX82ipnH8/U83s+y+1Rod+ExBICJpV1OWzwN/djn3vHcOW/Z3sORLv+OHa/f4XVZgKQhExBdmxm0t9fzqrquYW1/Gf3vkJT798AZO9fb7XVrgKAhExFcTywr43kcW8VfXT+PRF/dyy5efZeuBk36XFSgKAhHxXThk3PW26XzvjoWcON3Pu7/2e1ZtOuB3WYGhIBCRjPGWqZX8/FOLmVZdzMe+u46v/qZVzxykgYJARDLKhNJ8fvCxK1h66STuWbWVv/7BBnr6o36XNabppXMiknHyc8J86X2XMn1CCfes2kr7yR5W3t5CcZ7+ZHlBIwIRyUhmxp3XTuXfbpvL8zuPsnzlao509vhd1pikIBCRjHbr/DpW/ul8th08yW3feI69x0/7XdKYoyAQkYx3/Zsm8OBHFtLe2cPylavZpzAYVQoCEckKl08Zx4N3LORYVy/L71vNgRPdfpc0ZigIRCRrzK0v59t3LOBIZzwMDnUoDEaDgkBEssplDRV8+88v51BHN+//5vMcP9Xrd0lZz7MgMLN7zOwVM3vJzB41s/JBy+42s1Yz22pmN3pVg4iMTfMnj+ObH7qc3UdO8ZFvr6W7T88ZpMLLEcETwGzn3BxgG3A3gJnNBJYBs4AlwNfMTN9dJyJJueKS8XzxfZeybvcxPvX9F4nqew1GzLMgcM79yjk38BrB1UBdYnop8JBzrsc5txNoBRZ4VYeIjF3vmDORz908k19tPsjfP77J73KyVroe0/tz4AeJ6VriwTCgLdEmIpK0D7+1kf0nuln5zA6m15Twp4sm+11S1kkpCMzsSaBmmEWfdc49lljns0A/8L2BzYZZ/5wxnZmtAFYANDQ0pFKmiIxxf7ukmdZDnfzD45uYVl3MoqbxfpeUVVI6NeScu8E5N3uYz0AIfAi4GfiAe/0Vgm1A/aDd1AH7htn3Sudci3OupaqqKpUyRWSMC4eMLy27lMnjC/nL761nz9FTfpeUVby8a2gJ8LfAO51zg/9feRxYZmZ5ZtYITAPWeFWHiARDaX4O993eQl80xorvrtOdREnw8q6hrwAlwBNmtsHMvgHgnNsEPAxsBn4J3Omc0/9jIpKypqpi7l02jy37O/j8zzb7XU7W8OxisXNu6hss+wLwBa+OLSLBdW1zNR+7qol/f2YHV1wynpvnTPK7pIynJ4tFZMz5mxtncFlDOZ/50cvsOtLldzkZT0EgImNOTjjEvcvnEQ4Zn/iPF+mLxvwuKaMpCERkTKqrKOSfb53Dy3tP8JWnWv0uJ6MpCERkzFoyu4b3zKvlK79p5eW2E36Xk7EUBCIypv3dLbOoLM7lv/5wg24pPQ8FgYiMaWWFOfzzrXPYdrCTLz6xze9yMpKCQETGvGtmVLN8QQP3/W6HThENQ0EgIoFw9580M744j//x6Mt6ZfUQCgIRCYTS/Bw+d/NMXt57ggdX7/K7nIyiIBCRwLh5zkSunFbJPau2clDfd3yGgkBEAsPM+PzS2fRGY/zjT/UuogEKAhEJlCmVRXzi2qn8/OX9/GH7Yb/LyQgKAhEJnBVXNVFbXsDnf7ZFF45REIhIAOXnhPnMTc1s2d/BI+v2+F2O7xQEIhJIN8+ZyPzJFdyzahudPf1+l+MrBYGIBJKZ8b9unsnhzh6+/ttgv5ROQSAigXVpfTlLL53E/c/u5FCAbydVEIhIoN11w3T6o46v/ia4owIFgYgE2pTKIm5rqec/1uxmz9FTfpfjCwWBiATep66fiplx769f9bsUXygIRCTwJpYVcPuiyfxofRuthzr9LiftFAQiIsDHr7mEgpwwX34qeKMCBYGICDC+OI8PLprMT/+4j9cOd/ldTlopCEREEu64spFIOMQ3nt7udylppSAQEUmoLsln2eX1/Gh9G/uOn/a7nLRREIiIDLLiqiacg5XP7PC7lLRREIiIDFJXUci75tXy0Au7OdzZ43c5aaEgEBEZ4uPXXEJ3X4zvPheMr7RUEIiIDHFJVTHXN1fz4OpddPdF/S7HcwoCEZFh3LG4kSNdvTy2Ya/fpXhOQSAiMowrLhlPc00J9z+7E+fG9reYeR4EZvY3ZubMrDIxb2Z2r5m1mtlLZnaZ1zWIiCTLzPjIlU1sO9jJs61j+7uNPQ0CM6sH3gbsHtR8EzAt8VkBfN3LGkRERuqWuROpLM7j/md3+l2Kp7weEXwR+O/A4HHVUuA7Lm41UG5mEz2uQ0QkaXmRMLdfMZnfbm1n5xh+7YRnQWBm7wT2Ouf+OGRRLTD426LbEm0iIhln2eX1RELG99fsvvDKWSqlIDCzJ81s4zCfpcBngc8Nt9kwbedciTGzFWa21szWtre3p1KmiMiIVZfm8/ZZE/jh2j1j9lbSlILAOXeDc2720A+wA2gE/mhmrwF1wHozqyE+AqgftJs6YN8w+17pnGtxzrVUVVWlUqaISEo+sHAyx0718cuNB/wuxROenBpyzr3snKt2zk1xzk0h/sf/MufcAeBx4PbE3UOLgBPOuf1e1CEiMhquaBpPY2UR33t+bD5p7MdzBP9JfMTQCtwH/KUPNYiIXLRQyHj/ggZeeO0YWw+c9LucUZeWIEiMDA4npp1z7k7n3CXOuTc759amowYRkVTcOr+O3EhoTI4K9GSxiMhFGFeUy02za/jJi3vH3EVjBYGIyEW6bX49Hd39/HrLIb9LGVUKAhGRi3TFJeOZVJbPI+v2XHjlLKIgEBG5SOGQ8Z7L6nh6WzuHOrr9LmfUKAhERJJw6/w6Yg4efXHsvJ5aQSAikoTGyiLmT67gkXVtY+b11AoCEZEkvXd+Ha8e6uSlthN+lzIqFAQiIkl6x5yJ5EVC/Hh9m9+ljAoFgYhIkkrzc7iuuZqfv3yAaCz7Tw8pCERERuCdcydxuLOH1TuO+F1KyhQEIiIjcG1zNUW5YX76x3Nenpx1FAQiIiOQnxPm7bNq+MXGA/T2x/wuJyUKAhGREbpl7kROnO7jd69m95dnKQhEREZo8dQqygpysv70kIJARGSEciMhbppdwxObD3K6N3vfSKogEBFJwS1zJ9HVG+Xpbdn7RlIFgYhIChY2jqO8MIdVmw76XcqIKQhERFIQCYe4vnkCv95ykL5odt49pCAQEUnRktk1dHT3Z+3DZQoCEZEUXTmtksLcML/ceMDvUkZEQSAikqL8nDBXT6/iic0HiWXhu4cUBCIio+DGWTUcOtnDi3uO+11K0hQEIiKj4NrmaiIh41ebsu/0kIJARGQUlBXk8JaplazadCDrvrlMQSAiMkreNnMCrx05xfb2Lr9LSYqCQERklFzXXA3Ab7dm11PGCgIRkVFSW17AjAklPPWKgkBEJLCuaa5izc6jnOzu87uUi6YgEBEZRdfNqKY/5vh962G/S7loCgIRkVE0f3IFpfmRrDo95GkQmNknzWyrmW0ys38Z1H63mbUmlt3oZQ0iIukUCYe4anoVv9nanjVPGXsWBGZ2LbAUmOOcmwX8a6J9JrAMmAUsAb5mZmGv6hARSbfrmqtpP9nD5v0dfpdyUbwcEXwc+CfnXA+Ac25gnLQUeMg51+Oc2wm0Ags8rENEJK2unl6FGVlzesjLIJgOXGlmz5vZ02Z2eaK9FtgzaL22RJuIyJgwvjiPuXXlWfM8QSSVjc3sSaBmmEWfTey7AlgEXA48bGZNgA2z/jkn0sxsBbACoKGhIZUyRUTS7spplXztt9vp6O6jND/H73LeUEojAufcDc652cN8HiP+L/0fu7g1QAyoTLTXD9pNHbBvmH2vdM61OOdaqqqqUilTRCTtFk+tJBpzPLc987+sxstTQz8BrgMws+lALnAYeBxYZmZ5ZtYITAPWeFiHiEjazWuooDA3zLOvZv7zBCmdGrqAB4AHzGwj0At8yMVfybfJzB4GNgP9wJ3OuaiHdYiIpF1uJMTCxnFZ8WCZZ0HgnOsFPnieZV8AvuDVsUVEMsHiaVX8Zutm9h4/TW15gd/lnJeeLBYR8ciV0yoBePbVdp8reWMKAhERj0yrLqa6JI/fZfh1AgWBiIhHzIzFUyv5w/YjGf26CQWBiIiHFk+r5GhXb0a/bkJBICLiocVT49cJMvnuIQWBiIiHqkvzuaSqiNU7MvfBMgWBiIjHFjaNZ+1rx4hm6HUCBYGIiMcWNo7jZE8/m/dl5nUCBYGIiMcWNY0HyNjTQwoCERGPTSjNp7GyiOd3KghERAJrYeM41uw8mpHXCRQEIiJpsKhpPB3d/WzJwOcJFAQiImmwsGkcAM/vPOpzJedSEIiIpMHEsgImjy/k+Qy8YKwgEBFJk4WN41jz2tGMe++QgkBEJE0WNo7n+Kk+th486XcpZ1EQiIikyYLG+HWCdbuO+VzJ2RQEIiJpUldRQFVJHusVBCIiwWRmzG+oYK2CQEQkuOZPrmD30VMcOtntdylnKAhERNJo/pQKANbvOu5zJa9TEIiIpNGsSaXkRkKs25U5D5YpCERE0igvEmZObVlG3TmkIBARSbP5kyvYuLeD7r6o36UACgIRkbSbP7mC3miMTftO+F0KoCAQEUm7yybHLxivfS0zTg8pCERE0qyyOI8p4wsz5jqBgkBExAeXTa5g/e5jOOf/C+gUBCIiPpjXUMHhzl72Hj/tdykKAhERP8ytKwPgj3v8v2CsIBAR8UFzTSm54RAvtfn/hLFnQWBml5rZajPbYGZrzWxBot3M7F4zazWzl8zsMq9qEBHJVLmREG+aVMqGPWM4CIB/Af7BOXcp8LnEPMBNwLTEZwXwdQ9rEBHJWJfWlfHy3hNEff7GMi+DwAGliekyYF9ieinwHRe3Gig3s4ke1iEikpHm1pdzqjfK9vZOX+uIeLjvvwZWmdm/Eg+ctyTaa4E9g9ZrS7Tt97AWEZGMM6euHIANe44zfUKJb3WkNCIwsyfNbOMwn6XAx4G7nHP1wF3A/QObDbOrc8ZFZrYicW1hbXt7eypliohkpKbKIkryIr5fME5pROCcu+F8y8zsO8BfJWZ/CHwzMd0G1A9atY7XTxsN3vdKYCVAS0uL/09ciIiMslDIeHNdme+3kHp5jWAfcHVi+jrg1cT048DtibuHFgEnnHM6LSQigTS3vpxXDvj7JlIvrxF8FPh/ZhYBuonfIQTwn8CfAK3AKeDDHtYgIpLR5taV0xd1bNnfwbyGCl9q8CwInHPPAvOHaXfAnV4dV0Qkm8ytH3jC+LhvQaAni0VEfFRTmk91SR4vtfl3nUBBICLiIzPjzbVlbPTxS2oUBCIiPps1qZTWQ52c7vXngrGCQETEZ7Nqy4g5eOVAhy/HVxCIiPhs1qT423g27VMQiIgEUm15AWUFOb59mb2CQETEZ2bG7NpSjQhERIJs1qQyXtl/kr5oLO3HVhCIiGSAWZNK6Y3GaD2U/ldSKwhERDLArEnxJ4z9OD2kIBARyQCNlUUU5ITZuDf9F4wVBCIiGSAcMt40sYTNGhGIiATX7NoyNu07QSzN32GsIBARyRCzJpXS1Rtl19FTaT2ugkBEJEMMXDBO93UCBYGISIaYNqGYSMjS/s4hBYGISIbIi4Rpqirilf0n03pcBYGISAaZUVPKKwcUBCIigdVcU8Le46fp6O5L2zEVBCIiGaS5pgSAbWkcFSgIREQySPPE+HcTpPP0kIJARCSDTCrLpyQ/ktY7hxQEIiIZxMxorilhq0YEIiLBNaOmhFcOnMS59LxqQkEgIpJhZtSUcrK7n30nutNyPAWBiEiGeVPizqGtabpOoCAQEckw0xNBsCVNTxgrCEREMkxpfg615QVpu2CsIBARyUDNNSVpu4VUQSAikoGaJ5awo72L3v6Y58dKKQjM7DYz22RmMTNrGbLsbjNrNbOtZnbjoPYlibZWM/tMKscXERmrZtSU0h9zbG/v9PxYqY4INgLvAZ4Z3GhmM4FlwCxgCfA1MwubWRj4KnATMBNYnlhXREQGmTEh8c6hg95fJ4iksrFzbgvEn4QbYinwkHOuB9hpZq3AgsSyVufcjsR2DyXW3ZxKHSIiY01jZRGRkKUlCLy6RlAL7Bk035ZoO1+7iIgMkhsJMaWyiG0HvT81dMERgZk9CdQMs+izzrnHzrfZMG2O4YNn2GeozWwFsAKgoaHhQmWKiIw5c+vK6erp9/w4FwwC59wNI9hvG1A/aL4O2JeYPl/70OOuBFYCtLS0pOeFGyIiGeTf/svctBzHq1NDjwPLzCzPzBqBacAa4AVgmpk1mlku8QvKj3tUg4iIXISULhab2buBLwNVwM/NbINz7kbn3CYze5j4ReB+4E7nXDSxzSeAVUAYeMA5tyml/wUiIpISS9drTlPR0tLi1q5d63cZIiJZxczWOedaLrSeniwWEQk4BYGISMApCEREAk5BICIScAoCEZGAy4q7hsysHdg1qKkMOJHEfCVw2KPyhh5rtLa50DrnWz5cexD660Lrqb+SWy+V/hrapv5Kvm3wfCr9Ndk5V3XBtZxzWfcBViY5vzZdtYzWNhda53zLh2sPQn9daD31V/r6a2ib+iu13zkv+2vgk62nhn6a5LyXRnKsi9nmQuucb/lw7UHorwutp/5Kbr1U+mtom/or+bZ09ll2nBpKlZmtdRfxUIXEqb+So/5KjvorOenor2wdESRrpd8FZBn1V3LUX8lRfyXH8/4KxIhARETOLygjAhEROQ8FgYhIwCkIREQCLtBBYGbvMrP7zOwxM3u73/VkAzNrMrP7zewRv2vJRGZWZGbfTvxefcDverKBfqeS48XfrawNAjN7wMwOmdnGIe1LzGyrmbWa2WfeaB/OuZ845z4K/BnwPg/LzQij1Gc7nHN3eFtpZkmy394DPJL4vXpn2ovNEMn0WRB/p4ZKsr9G/e9W1gYB8C1gyeAGMwsDXwVuAmYCy81sppm92cx+NuRTPWjT/5nYbqz7FqPXZ0HyLS6y34h/D/eexGrRNNaYab7FxfeZjKy/Ru3vVkpfVekn59wzZjZlSPMCoNU5twPAzB4Cljrn/g9w89B9mJkB/wT8wjm33tuK/TcafRZEyfQb0EY8DDaQ3f/QSkmSfbY5vdVlnmT6y8y2MMp/t8baL2otr/9rDOL/Uda+wfqfBG4A3mtmf+FlYRksqT4zs/Fm9g1gnpnd7XVxGex8/fZj4FYz+zppfk1AFhi2z/Q7dV7n+x0b9b9bWTsiOA8bpu28T8w55+4F7vWunKyQbJ8dAYIamoMN22/OuS7gw+kuJkucr8/0OzW88/XXqP/dGmsjgjagftB8HbDPp1qyhfpsZNRvyVOfJSdt/TXWguAFYJqZNZpZLrAMeNznmjKd+mxk1G/JU58lJ239lbVBYGbfB54DZphZm5nd4ZzrBz4BrAK2AA875zb5WWcmUZ+NjPoteeqz5PjdX3rpnIhIwGXtiEBEREaHgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAJOQSAiEnAKAhGRgFMQiIgE3P8Hn0l220nldZ8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogx(omega, numpy.unwrap(numpy.angle(G)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}